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We study the patterns formed by adding N sand-grains at a single site on an initial periodic 
background in the Abelian sandpile models, and relaxing the configuration. When the heights at all 
sites in the initial background are low enough, one gets patterns showing proportionate growth, with 
the diameter of the pattern formed growing as N l ^ d for large N, in d-dimensions. On the other hand, 
if sites with maximum stable height in the starting configuration form an infinite cluster, we get 
avalanches that do not stop. In this paper, we describe our unexpected finding of an interesting class 
of backgrounds in two dimensions, that show an intermediate behavior: For any N, the avalanches 
are finite, but the diameter of the pattern increases as N a , for large N, with 1/2 < a < 1. Different 
values of a can be realized on different backgrounds, and the patterns still show proportionate 
growth. The non-compact nature of growth simplifies their analysis significantly. We characterize 
the asymptotic pattern exactly for one illustrative example with a = 1. 
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I. INTRODUCTION 

In the last two decades a large amount of study has 
been devoted to understanding various models of self- 
organized criticality, in particular, the Abelian Sandpile 
Model (ASM) (see [IJ[2] f° r reviews). These have mainly 
dealt with the critical exponents of avalanches produced 
in sandpiles driven slowly in their critical steady state. 
But the ASM has other interesting properties, not di- 
rectly related to its critical exponents. In particular, one 
sees very interesting and beautiful spatial patterns when 
many sand grains are added at a single point on an ini- 
tially periodic background, and we relax the configura- 
tion using the ASM toppling rules. One such pattern on 
a square lattice is shown in Fig. [T] 

Our interest in these patterns comes from several rea- 
sons. Firstly, these are analytically tractable examples 
of complex patterns that are obtained from simple de- 
terministic evolution rules. Of course, there are many 
known examples of complex patterns obtained this way 
(e.g. Conway's game of life [3]). But a detailed character- 
ization of such patterns is usually not easy. The sandpile 
patterns studied here are special, as they are nontriv- 
ial, but of intermediate complexity, and are analytically 
tractable. 

Secondly, growing sandpiles studied here are qualita- 
tively different from the growth models that have been 
studied in physics literature earlier, such as the Eden 
model, the diffusion limited aggregation, or the sur- 
face deposition (3HS]. These are the simplest models of 
proportionate growth, a well-known feature of biological 
growth in animals, where different parts of a growing an- 
imal grow at roughly the same rate, keeping their shape 
almost the same. In the models of growth studied earlier, 



growth is confined to some active outer region. The in- 
ner structures, once formed are frozen, and do not evolve 
further in time. This is not the case for the patterns 
studied in this paper. Figure [I] shows two patterns pro- 
duced on the same background but with different values 
of the number A" of grains added. As A" increases, the 
pattern grows in size, but we see that while some new de- 
tails emerge near the center, the relative proportions of 
different parts in the outer region of the pattern remains 
unchanged. As the pattern grows, different features, not 
only grow in size, but also are moved in space with time. 

The third motivation for our study is some intriguing 
connection of these to the mathematics of discrete ana- 
lytic functions. We have not explored this much, but will 
discuss briefly in an appendix. 

There have been several earlier studies of the spatial 
patterns in sandpile models. First of them was by Liu 
et.al. [7]. The asymptotic shape of the boundaries of 
the patterns produced in centrally seeded sandpile model 
on different periodic backgrounds was discussed in 
Borgne et. al. |9 obtained bounds on the rate of growth of 
these boundaries, and later these bounds were improved 
by Fey et.al. [10] and Levine et.al. [11 . A detailed 
analysis of different periodic structures found in the pat- 
terns were first carried out by Ostojic [12] who also first 
noted the exact quadratic nature of the toppling function 
within a patch. Wilson et.al. [13] have developed a very 
efficient algorithm to generate patterns for a large num- 
bers of particles added, which allows them to generate 
pictures of patterns with A" up to 2 26 . 

Other special configurations in the Abelian sandpile 
models, like the identity [9j [HI [15] or the stable state 
produced from special unstable states, also show complex 
internal self-similar structures [7], which share common 
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FIG. 1. (Color online.) The stable configurations for the 
Abelian sandpile model on a square lattice, obtained by 
adding N grains at the origin. In the initial configuration 
all heights are 2. For comparison the size of the first figure 
has been enlarged by a factor 2. (Details can be seen in the 
electronic version using zoom in). 



features with the patterns studied here. In particular, 
the identity configuration on the F-lattice has recently 
been shown to have similar spatial structures [T5] . 

There are other models, which are related to the 
Abelian sandpile model, e.g., the Internal Diffusion- 
Limited Aggregation (IDLA) 16 , Eulerian walkers (also 



called the rotor-router model) |17IfT9] . and the infinitely- 
divisible sandpile [11], which also show similar structure. 
For the IDLA, Gravner and Quastel showed that the 
asymptotic shape of the growth pattern is related to the 
classical Stefan problem in hydrodynamics, and deter- 
mined the exact radius of the pattern with a single point 
source [20] . Levine and Peres have studied patterns with 
multiple sources in these models, and proved the exis- 
tence of a limit shape[21j. Limiting shapes for the non- 
Abelian sandpile has recently been studied by Fey. et.al. 
[22]. 

The standard square lattice produces complicated pat- 
terns and it has not been possible to characterize them 
fully, so far. In an earlier paper [23], we considered the 
pattern produced on a F-lattice (Fig. [2^a)), and deter- 
mined exactly the sizes of different patches in the asymp- 
totic pattern. The pattern produced by adding grains at 
one site on a background with a periodic chequerboard 
pattern of alternate sites with heights 1 and 0, is shown 
in the Fig. [2|b). In paper [23], we studied the patterns 
when sink sites are present, or when addition is made at 
more than one site. In paper 25\, we have studied the 
effect of noise on such patterns. 

If the average initial height in a background is high, one 
gets infinite avalanches, with the diameter of the pattern 
becoming infinite for finite number of particles added. 
Such backgrounds have been termed as 'explosive'. In 
other cases, the diameter of the pattern is finite for any 
finite N, and increases as N 1 ^ in d- dimensions. We call 
such a growth as compact growth. All the patterns stud- 
ied in [23H25] showed compact proportionate growth. In 
this paper, we describe a remarkable class of patterns 
where the diameter remains finite for any finite N, but 
grows as A^ a , with 1/2 < a < 1. We call such growth as 
non-compact proportionate growth. Characterization of 
these patterns, as will be shown, is simpler than in the 
compact case. 

We found two classes of backgrounds, both infinite, on 
a directed triangular lattice (see Fig. EL for which the 
growth is proportionate, with the growth exponent a — 
1. Our numerical study shows, but we have no formal 
proof, that different backgrounds belonging to the same 
class produce the same asymptotic pattern. In addition, 
we found infinitely many backgrounds on the F-lattice 
which produce patterns with proportionate non-compact 
growth. However, in these cases the growth exponent 
a takes a different value, with 1/2 < a < 1 for each 
member. 

There are some earlier works on the growth rate of the 
sandpile patterns. Different non-explosive backgrounds 
for a DASM were studied in [TU] Q~T] [2E] . However, in all 
these examples, studied so far, the growth of the patterns 
is compact. For a DASM on a square lattice, it was 
shown [26], that the pattern produced on a background 
of constant height z < z c — 2, is always enclosed inside a 
square whose width grows as yN. 

We also discuss the exact characterization of the pat- 



tern shown in Fig. 11 one of the two asymptotic patterns 
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FIG. 2. (Color online.) (a) A F-lattice with directed bonds, 
and a chequerboard distribution of grains on it. Unfilled cir- 
cles denote height z = I and filled ones z = 0. (b) The stable 
configuration obtained by adding 10 5 grains at the origin. 
Color code: red =0, white = 1. The apparent orange regions 
in the picture represent the patches with checkerboard config- 
uration. (Details can be seen in the electronic version using 
zoom in.) 



we have found with a — 1. This is described, as in the 
earlier studied case of compact growth [23], in terms of 
the scaled toppling function. However, the analysis of 
non-compact patterns is actually simpler. Clearly, for 
a > 1/d, the mean excess density of particles in the 
toppled region is zero, for the asymptotic non-compact 
growth patterns. The patterns are made of large re- 



gions where the heights are periodic, and we call them as 
patches. We find that, inside each patch, the mean den- 
sity is exactly the same as in the background, and the ex- 
cess grains are concentrated along the patch boundaries. 
There are also some boundaries where excess grains den- 
sity is negative. We show that this leads to the scaled 
toppling function being a piece- wise linear function of the 
rescaled coordinates. Thus, in each patch, the potential 
function is specified by only 3 coefficients. In contrast, 
for the compact patterns, the scaled toppling function is 
a quadratic function of the coordinates in each patch [23] , 
and one has to determine 6 coefficients for each patch, to 
determine the function fully. 

We are able to reduce the problem of determining the 
asymptotic pattern in Fig. 11 to that of finding the lat- 
tice Green's function on a hexagonal lattice. This is 
known to be expressible as integrals that can be eval- 
uated in closed form [27], and this leads to a full solution 
of the problem. For the compact growth patterns studied 
earlier on F-lattice (23] [2J] , we define a discrete analytic 
function F p (z), which is the discrete analogue of the an- 
alytic function z p , for any rational value of p, and show 
that the patterns can be characterized in terms of this 
function. 

The plan of this paper is as follows: In section [IT] we 
discuss how different periodic background configurations 



give rise to different rates of growth. In section[TTT] we de- 
scribe two classes of periodic background configurations 
on a directed triangular lattice that give rise to patterns 
which show proportionate growth with a — 1. In sec- 



tion IV we argue that, for any pattern with non-compact 
proportionate growth the rescaled toppling function is 
piece- wise linear. In section [V] we discuss exact char- 
acterization of the simplest of the patterns with a = 1. 
Patterns on the F-lattice, with 1/2 < a < 1 are discussed 
in section |VI[ Section |VII| contains some concluding re- 
marks and a discussion about connection to tropical poly- 
nomials. An appendix describes how the characterization 
of these patterns involves the theory of discrete analytic 
functions defined on many-sheeted discretized Riemann 
surfaces. 



II. THE COMPACT AND NON-COMPACT 
GROWTH PATTERNS 

The simplest growing patterns are found in the Manna- 
type sandpile models with stochastic toppling rules [28] . 
In these models, when the density of particles p Q in the 
background is small, the avalanches are always finite. In 
the relaxed configuration, the toppled sites form a nearly 
circular region (see Fig. pjj). The asymptotic pattern 
seems to be perfectly circular disc of uniform density, 
with an average density p* inside the circle and p Q out- 
side. The value of p* is independent of the background 
density p a , and is equal to the unique steady state density 
of the corresponding self organized critical model with 
random sites of addition, and dissipation at the bound- 



FIG. 3. The pattern produced by adding N = 10 5 grains at a 
single site on a stochastic ASM defined on an infinite square 
lattice, and relaxing the configuration. The initial configura- 
tion has all sites empty. The threshold height z c = 2, and 
on toppling two grains are transfered either to the vertical 
or horizontal neighbors, with equal probability. Color code: 
White=0, and Black=l. 



ary [28]. The region inside the circle forgets about the 
initial height configuration, and is in the self-organized 
critical state. The boundary of the affected region is thin 
with a sharp transition of density from p* to p G . Then 
considering that, for large N, all the added grains are 
confined inside the circular region of diameter 2A, we get 

N = (p* — p ) 7rA 2 + Lower order in A. (1) 

Thus the pattern has a compact growth. 

For densities p Q close to, but below p*, sometimes a 
single particle addition can lead to very large increase in 
the size of the toppled region. However, probability of 
such large jumps decreases exponentially with size, and 
for any finite N, with probability 1, avalanches remain 
finite. As long as p a is less than p*, the system relaxes, 
forming a pattern whose diameter grows as y/N. Adding 
a single grain on a background of super-critical density 
(p > p*) gives rise to infinite avalanches, with non-zero 
probability. Then, with probability 1, such backgrounds 
will lead to an infinite avalanche for some finite value 
of N. In higher dimensions also, a similar behavior is 
expected. 

In the models with deterministic relaxation rules there 
is no simple quantifier like critical density p*, separating 
the explosive and non-explosive backgrounds. Whether 
the periodic or random background is explosive or not de- 
pends in a complicated way on the built-in height correla- 
tions. For example consider the BTW model on a square 
lattice, where the steady state density p* = 17/8 = 2.125 
[28] . It has been shown that a background with a ran- 
dom assignment of height 3 with probability e, on a sea 



FIG. 4. A directed triangular lattice. 



of constant height 2 is explosive, even for arbitrary small 
value of e [26] , although the average density p Q — 2 + e is 
much less than p*. On the other hand, it is also possible 
to construct a robust periodic background with density 
arbitrarily close to maximum stable height 3 [26] . 

We will show that there is a large class of backgrounds, 
with a range of densities, for which the growth is less than 
explosive, but more than compact. 



III. EXAMPLES OF NON-COMPACT GROWTH 

We first discuss the patterns with a — 1. We start with 
a ASM on a directed graph corresponding to an infinite 
two dimensional triangular lattice, with each site having 
three incoming and three outgoing arrows (see Fig. [4|. 
The threshold height z c — 3, for each site. If the height 
at any site is above or equal to z c , it is unstable, and 
relaxes by toppling: in each toppling, three sand grains 
leave the unstable site, and are transferred one each along 
the directed bonds going out of the site. 

We consider two classes of backgrounds on this lattice: 

Class I: We consider the lattice as made of trian- 
gular plaquettes, which are joined together to make 
tiles in the shape of regular hexagons with edges of 
length We cover the two-dimensional plane with 
these tiles. Sites that lie on the boundaries of these 
hexagons are assigned height 1, and the rest of the 
sites have height 2. Figure [5ja) shows the back- 
ground configuration for the case 1 — 2. 

Class II: For these backgrounds, we cover the two- 
dimensional plane with tiles in the shape of equilat- 
eral triangles of edge- length The sites that lie on 
the boundaries of the triangles, and are shared by 





FIG. 6. (Color online.) The class I pattern, formed by adding 
N = 500 particles at the origin on the first background shown 
in Fig. [5} The apparent uniform green color of the back- 
ground is actually a periodic structure. Details can be viewed 
in the electronic version using zoom in. 
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FIG. 5. Examples of backgrounds of class I and II, respec- 
tively. The filled circles represent height 1 and unfilled ones 
height 2. 



two triangles, are assigned height 1, and remain- 
ing sites are assigned height 2. Sites that are at 
the corners of triangular tiles, and shared by six of 
them, are also assigned heights 2. The background 
configuration corresponding to I — 4 is shown in 
figure [5^6). The pattern made of triangular tiles 
with / = 3 is same as the class I background with 
hexagon of edge- length 1. Hence, only patterns 
formed with triangles of edge-length / > 4 will be 
said to be in this class. 

The patterns produced by adding N grains, where N 
is large, at a single site on the two backgrounds in Fig. 
[5] are shown in Fig. [6] and [7] While the patterns look 
quite similar, a closer examination shows that they are 
not identical. In Fig. [7| there are extra lines of particles 



FIG. 7. (Color online.) The class II pattern, formed by adding 
N = 500 particles at the origin on the second background in 
Fig. [5] The apparent uniform green color of background is 
actually a periodic structure. Details can be viewed in the 
electronic version using zoom in. 
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FIG. 8. The patterns in terms of Q(r), corresponding to 
those in Fig. [6] and [7] Sites with zero Q (r) are colored white, 
and non-zero are colored black. The larger patches are given 
identifying labels. 



within the brownish patches which break each patch into 
smaller parts. In fact, with the identification of some 
patches having only a point in common, as discussed 
later, we can show that each patch breaks into exactly 
three patches. These three parts have similar periodic 
pattern, but with different orientations. 




2000 4000 6000 8000 
N 

FIG. 9. (Color online.) The diameter 2A of the patterns as a 
function of the number N of added grains. The cases shown 
are (i) Class I, I = 1, (ii) Class I, I = 2, and (hi) Class II , 
1 = 4. The corresponding straight line fits have slopes given 
by 1.1, 2.7 and 1.7 respectively. 



This differences can be seen more clearly in terms of the 
net excess change in height Q (R) in a unit cell centered at 
R, where the unit cell is that of the background pattern. 



Q(R)= J2 + (2) 

i?'Gunit cell 



where Az (R) is the change in height at site R. For ex- 
ample in the first background in Fig. [5] a unit cell is a 
hexagon of edge length 1 = 2, and for the second back- 
ground it is a parallelogram of each side length I = 4. A 
site that is on the edge of the unit cell is counted with 
weight 1/2, and a site on the corner of the hexagon with 
weight 1/3, and on the corner of the parallelogram with 
weight 1/4. By construction, the function Q (r) is zero 
inside each patch, and non-zero along the boundaries be- 
tween patches. The patterns in terms of these variables, 
corresponding to those in Fig. [6] and [7] are shown in Fig. 

El 

We have seen that the patterns on these two classes of 
backgrounds exhibit proportionate growth, i.e., all the 
spatial features inside the patterns for large N, grow at 
the same rate with the diameter. We define the diameter 
2 A, in general, for any pattern in this paper, as the height 
of the smallest rectangle containing it. For the patterns 
in Fig. [6] and Fig. [7| it is then the length of a side of 
the bounding equilateral triangle. This particular choice 
makes 2A as an integer multiple of \/3, on the triangular 
lattice. We find that for both types of backgrounds in 
Fig. [5} the diameter of the pattern grows linearly with 
N (Fig. g. 
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IV. PIECE- WISE LINEARITY OF THE 
TOPPLING FUNCTION. 

Considering the proportionate growth, let us define a 
rescaled coordinate r — R/N a , where R = (x,y) is the 
position vector of a site on the lattice. The number of 
topplings at any site inside the pattern, scales linearly 
with N. Let us define 



(r) — lim — — — - 

V ' N^oo N 



(3) 



We now show, using an extension of the argument given 
in [23], that the function <fi is linear inside periodic 
patches in all the patterns with non-compact growth, 
i.e., with a > 1/2. Within a patch, the function 0(f) 
is expandable in Taylor series around any point f , not 
on the boundary of the patch. Defining f a = (£ G , r) ), and 
Af G = (A£, Ar)) we have 



(£,+ A£, r7 o + Ar7)-0(£ o ,r7 o ) 

= dA£ + eAr/ + (9 (A£ 2 , Ar/ 2 , A^Ary) 



(4) 



Consider any term of order > 2 in the expansion, for 
example, the term ~ (A£) 2 . This can only arise due to 
a term ~ (Ax) 2 7V 1_2a in the toppling function Tn(R). 
Then, considering the fact that Tn(R) is an integer func- 
tion of x and y, it is easy to see that this term would 
lead to discontinuous changes in Tn(R) at intervals of 
Ax ~ 0(N a ~ 1 / 2 ). As a > 1/2 for non-compact growth 
patterns, this leads to a change in the periodicity of 
heights at such intervals inside each patch which them- 
selves are of size ~ N a . This would then result in many 
defect lines within a patch, in the pattern at large N. 



However there are no such features in Fig. 11 Therefore 
inside each periodic patch, 0(f) must be exactly linear in 
f. In fact, it turns out that the integer toppling function 
Tjv(-R) is exactly linear inside a patch even for any finite 
N, except for an additional periodic term of periodicity 
equal to that of the heights inside the patch. 

Another consequence of the exact linearity of the po- 
tential function in each patch is that all patch boundaries 
in the asymptotic pattern are straight lines. 

The argument finally relies on the two observed (not 
rigorously established) features of the patterns, i.e., there 
is proportionate growth, and that the patterns can be de- 
composed in terms of periodic patches which are them- 
selves of size O (N a ). 

Let us write the toppling function Tjv (R) within a sin- 
gle patch P as 



Tn(R) — Ap + Kp • R + T per i dic(R) , 



(5) 



where T per i dic(R) is a periodic function of its argument 
with zero mean value. If e\ and e 2 are the basis vectors 
at the unit cell of the periodic pattern then we have 



T N {R + ei ) 
T N (R + e 2 ) 



T N (R)=K P -e 1 , 
T N {R)=K P -e 2 . 



As Tjy(R) are integer valued functions, Kp ■ e± and 
Kp ■ £2 can only take integer values. If g\ and g 2 are the 
unit vectors in the reciprocal space of the super lattice of 
the periodic pattern, 



9i 



hJ 



(7) 



then Kp must be an integer linear combination of g\ and 
g 2 , and can be written as 

Kp = mgi + n 2 <?2, (8) 
where n\ and are some integers. For example, in the 



background pattern in Fig. [TOj a choice of the basis vec- 
tors and its reciprocal vectors is 



9i 




e 2 



92 




(9) 



The fact that K p is constant inside a patch, implies that 
the patches can be labeled by the pair of integers (711,712). 

An interesting consequence of this linear dependence 
of Tn(R) is that there are no transient structures within 
the patches. On increasing N, if the A p function in- 
creases, all sites in the patch P, except possibly those 
at the patch boundaries, undergo same number of addi- 
tional topplings. 



V. CHARACTERIZING THE CLASS I 
ASYMPTOTIC PATTERNS 

We now discuss characterization of the asymptotic pat- 
tern of class I, showing a = 1. In this section we quanti- 
tatively characterize the asymptotic pattern for the case 
1 = 1. The background configuration is shown in Fig. 10 
A site on the triangular lattice can be labeled uniquely 
by a pair of integers (p,q), such that its position on a 
complex plane can be written as R = p + qoj, where 
uj = exp(z27r/3) is a complex cube root of unity. Then, 
the height variables in the background pattern in Fig. 
lOl can be written as 



h(p + quo) = 2ifp + g = (mod 3), 

= 1 otherwise. (10) 

The average height in the background, (z) =4/3. The 
configuration of the pile produced on this background, 
by adding TV = 3760 grains at the origin is shown in Fig. 

HU 

We see that the sites toppled due to addition of the 
grains are confined within an equilateral triangle. The 
pattern can be thought of as a union of patches, inside 
which the heights are periodic. A zoom-in showing the 
height configuration with five patches meeting at a point 



(6) is shown in Fig. 12 There are only two types of periodic 
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FIG. 12. (Color online.) An example of patch boundaries 
in Fig. 1 1 1 1 meeting each other. Each filled hexagon repre- 
sents Wigner cell around a site, and the color in them denotes 
height of that site. The color code is same as in Fig. |11| 



FIG. 10. The background of class I corresponding to Z = 1. 
The filled circles represent height z = 1 and unfilled ones 
z = 2. 
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FIG. 11. (Color online.) The pattern produced by adding 
N = 3760 grains at a single site on the background in Fig. 
|10| and relaxing. Details can be seen in the online version 
using zoom-in 



patches seen: one is like the background, where the sites 
of height 2 are surrounded by sites of height 1, and the 
other with heights surrounded by heights 2. Then, the 
average height (z) inside both types of patches are same. 
In fact, it is equal to that of the background, (z) — 4/3. 

The patches in the outer region of the pattern are 
big, and they become smaller, and more numerous as 



we go inwards. Along the common boundary of adja- 
cent patches, we see line-like defect structures, and only 
along these lines the density is different from the back- 
ground. In Fig. [12] one can also see the periodicity of 
the structures along the patch boundaries. Some patch 



boundaries, like the horizontal boundary in Fig. 12 have 



a deficit of particles compared to the background. 

The boundaries of the patches are seen more clearly in 
terms of Q(R) variables, as shown in figure [8ja), where 
we have labelled different patches as A, A',B,B'... etc.. 

The dependence of 2A on A" for this background is 
shown in Fig. [9} We see that the diameter for the pattern 
grows asymptotically linearly with N, but it grows in 
bursts: it remains constant for a long interval as more 
and more grains are added, and suddenly increases by a 
large amount at certain values of N. For example, at A" = 
3721, the 2A is 2276\/3, and it jumps to a value 2408\/3 
when one more grain is added. Let Jmax(N m ) denote 
the size of the maximum jump in 2 A encountered, as N 
is varied from 1 to N m . In Fig. 13 we have plotted the 
variation of Jmax(N m ) with N m . The graph is consistent 
with a power-law growth, with a power around 1/2. Thus 
the fractional size of the bursts decreases for large A". 

We define scaled complex coordinates r = H/N, where 
R = p + quj is the complex coordinate of the site (p, q). 
We define the rescaled toppling function for this pattern 
as 



(r) = lim 



V3T N (rN) 



2N 



(11) 



Then it is easy to see that = (d^4>, d v (j)) is equal to 
the mean flux of particles at r. If we consider a small 
line element dl = (d£,dr]), then the net flux of particles 
across the line dl equals NV(f)-dl. Then, the conservation 
of sand grains implies that the toppling function Tn (R) 
satisfies the equation 

V 2 Q T N (R) = Sz (R) - N5 (R) , (12) 

where V 2 , is the finite-difference operator on the lattice, 
corresponding to the Laplacian V 2 . It is easy to see that 
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FIG. 13. Jmax{N) for the background in Fig. 10 has a square 



root dependence on N with a fitting function 3.25-y/x. 




FIG. 14. The 1/R transformation of the pattern in Fig. 11 
where R is the complex conjugate of R. Labels are the same 
as used in Fig. |8[a). 



this implies that the scaled potential function <fi satisfies 
the Poisson equation 

V 2 0(r) = Ap(r) -6(r), (13) 

where Ap (r) is the areal density of excess grains at r. It 
is related to (Az (r)), the mean excess grain density per 
site by 



Ap(r) = -^(A«(r)). 



(14) 



The piece-wise linearity of (f) simplifies the analysis of 
the pattern, significantly. The potential function can be 
characterized by only three parameters. Using Eq. d8|, 
([9| and (11 ), for each patch P, we can find a pair of in- 




FIG. 15. The adjacency graph of patches in the pattern in Fig. 
1 1 1 1 The vertices corresponding to the brownish and greenish 
patches in the pattern are denoted by different colors. The 
pair of patches labeled by the alphabets and its corresponding 
primed alphabets in Fig. |8ja) are represented by same vertex 
on the graph. 



tegers (m, n) such that the potential in patch P is char- 
acterized by 



0(r) 



1 



2^3 



(■D m n r -|- ~D m n v^) -\- fm,7 



where 



D, 



(15) 



(16) 



and f m ,n is a real number, constant everywhere inside 
the patch. Here z denotes the complex conjugate of z. 

Each patch is characterized by a complex number D m n 
which is the coefficient in the potential function (f) (r) 
of the patch. In the complex D-plane, each patch with 
labels as in Fig. [8ja) can then be represented by a point. 
We connect two patches by a line if they share a common 
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boundary. Then the resulting figure, shown in Fig 
is the adjacency graph of the patches. 

We can determine the connectivity structure of this 
graph, without knowing the full potential function in 
each patch. We first take 1/R transformation of the 
This is shown in Fig 



pattern. 



14 



Some of the big- 
ger patches are denoted by capital alphabets in Fig. [8^a) 
and their corresponding patches on the transformed pat- 
tern in Fig. 14 The patches A and A' in Fig. 



|8j^a) are 

adjacent to the outer region O through the same verti- 
cal boundary. Matching the values of the function 0(r) 



10 



and fixing the discontinuity in its normal derivatives at 
the boundary, it is easy to see that 0(r) has the same 
functional form in the patches A and A'. In fact, it is 
convenient to imagine that the boundary between O and 
A moved to the right by an infinitesimal amount, so that 
it does not touch the patches D and I', and then A and 
A' would actually join to form a single connected patch 
A. We thus consider A and A' as one patch, and both 
can be represented as one point on the D-plane. Sim- 
ilarly, we identify B and B', C and C', etc. Then the 
adjacency graph can be constructed by joining the sites 
on the D-plane, according to the adjacency of patches in 
Fig. 
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It turns out that the patches corresponding to ra + n = 
2 (mod 3) do not appear in the pattern, and the adja- 
cency graph, as shown in Fig. |15[ is a hexagonal lattice 
with some extra edges shown in brown color. These extra 
edges connect all the vertices at same distance from the 
origin (0, 0) (in the L 1 metric), and also connect some of 
the diagonally opposite sites on the rectangular faces of 
the graph as shown in figure. 

The charge density Ap (r) is zero inside the patches, 
and the excess grains due to addition are distributed 
along the patch boundaries, leading to nonzero line 
charge densities separating neighboring patches. Then 
the density function Ap (r) is a superposition of the line 
charge densities along the patch boundaries. There are 
three kinds of line charges of charge density A = — l/y3, 
1, and 2/ VS. 

From the electrostatic analogy, it is seen that <fi (r) is 
continuous across the common boundary between neigh- 
boring patches, and its normal derivative is discontinuous 
by an amount equal to the line charge density A along the 
boundary. Let P and P' be the two neighboring patches 
with the equation of the boundary between them 



|r| exp (id) + A, 



(17) 



such that the patch P' is on the left of the boundary. 
Then using the continuity condition, it is easy to show 
that 



D, 



D, 



exp (id) and 



where A is the complex conjugate of A. We note that, 
there are only six different types of patch boundaries in 
the pattern, with angle 9 an integer multiple of 7r/6. 

It is easy to check that the matching conditions along 
the edges of hexagonal lattice (denoted by blue solid line 
in Fig. 15) are sufficient to determine -D m , n for all the 
vertices. The line charge density A = — l/v3 for the 
patch boundaries corresponding to these edges. Also, the 
potential function <fi — 0, for the vertex at the origin, and 
hence, D and /both vanishes. Then using the matching 
condition, it is easy to check that, the values of D m>n are 
consistent with the form in Eq. (16). 



The function / mjn satisfies the discrete Laplace's equa- 
tion on the underlying hexagonal lattice of the adjacency 
graph i.e. 



fm',n' - 3/ m>n = for (ra, n) ^ 0, (19) 



where (m',n f ) denotes the three neighbors of the vertex 
(m, n) on the hexagonal lattice. This can be checked 
from the concurrency condition of patch boundaries. For 
example consider the edges OA, DA' and FA on the 
adjacency graph. The corresponding patch boundaries 
in the pattern intersect at the same point (Fig. [8^ a)). 
Then it is easy to check using the matching condition in 



Eq.(18) that 



f p ,-f p = Re[A(p p , -DJ]/V3, 



(18) 



/o + /d + /i = 3/ a . (20) 



Similar equations hold for the other vertices. 

In the region outside the pattern, where none of the 
sites toppled, the potential function (f) (z) — 0. This cor- 
responds to m — n — 0, and /o,o = 0. The solution of 
the Laplace's equation with the above boundary condi- 
tion can be written in the following integral form [27] 



fm,n — 



47T 2 



1 — cos (ki (2m — n) /3 + A^n) 
1 — (cos 2&2 + 2 cos k\ cos fo) /3 



for m + n — (mod 3), 



(21) 



where / is a normalizing constant, which determines the function in region A, and C' is 
pattern up to a scale factor. For the sites with ra + n = 1 
(mod 3), f m , n are the average of those corresponding 
to the neighboring sites. As an example the potential 




(22) 
(23) 
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where r = £ + irj, and i — \J— 1. Then the equation of 
the patch boundary between patches A and O is 

t = //V3, (24) 

and that of the boundary between patches C' and O is 

V3£ + 3// + 2/ = 0. (25) 

Equivalently, the length of an edge of the bounding equi- 
lateral triangle of the pattern is equal to 2/A, for large 
A. 

The constant / in Eq. (21) can be calculated using 



1000 



the form of the potential function near the site of addi- 
tion. As noted, the function <fi can be considered as the 
potential due to line charges along the patch boundaries 
and a point charge of unit amount at the origin. Then, 
close to the origin the solution diverges logarithmically 
as (f) ( r ) — — (27r) — 1 log (|r|), and the potential function 
is an approximation to this solution by a piece-wise lin- 
ear function. Then, there are coordinates r D inside each 
patch (m,n) with \m\ + \n\ large, where the and its 
first derivatives are equal to and its first derivatives, 
respectively. Then, 



2^3 



2>/3^<t> (r)| ro 
or 

\^-^m,n^o "F I-^m,n^oj ~F fm,n 



1 

'2k 



The above two equations imply 



Sr. 



1 

2^ 



log (| 



hn and 

log (|r |) . 

(26) 

(27) 



for \m\ + \n\ large. Comparing it with the Eq. ( plj ) for 
large \m\ + |n| we find that the numerical constant / 
1 / \/3. This determines the potential function completely, 
and thus characterizes the pattern. For example, as in 
figure [8| a), the equation of the rightmost boundary of 
the pattern, using Eq. (24) is x — A/3. Equations of 



other boundaries of patches can be calculated similarly. 
For example, the reduced coordinates of the point where 
the patches D and D' meet in Fig. [8^a), is determined 
by the condition that it is a common point of patches _D, 
J and A' and that the function d> is continuous. 



1 



r i,o 



£ = f2,l 



— /3,1 




(28) 



Then using the values /i 5 o = 1/3-^, /2,i — l/2\/3, and 
/3 ( i = 7/6\/3 — 1/tt [27] we get the reduced coordinates 
of this point as 



fori) 



i 



i 



(29) 
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FIG. 16. The discrepancy 2AA between the actual height of 
bounding triangle, and the asymptotic value 2N/ a/3 plotted 
as a function of N. The straight line shows a simple power-law 
fit with power 3/4. 



Equivalently, the height of the bounding equilateral 
triangle increases as 2A/V3 ~ 1.154A. The estimated 
slope of the fitting line in Fig. [9] is 1.1, in reasonable 
agreement with the theory. However, even though the 
exact function A(A) has large fluctuations of number 
theoretic origin, the estimated slope is noticeably lower 
than the calculated asymptotic value. To examine this 
discrepancy closer, we have plotted in Fig. 16 the dis- 
crepancy 2AA = 2A/V3 - 2A(A) as a function of A. 
We find that this appears to increase with N as A 3 / 4 , 
for large A. The reason for this behavior is not under- 
stood yet. 

For backgrounds, with / > 1, our numerical results sug- 
gest that there is a crossover length R*(l), and initially, 
for R < R*(l), the avalanches grow "explosively" in size. 
As a result, the number of particles inside a disc of ra- 
dius R* in the final pattern is less than that in the initial 
background. The net flux of particles going out of the 
disc increases with R until the radius becomes of order 
R*. After this, the large-scale properties of the pattern 
are the same as that of I — 1 pattern, with the num- 
ber of particles added AiN, where Ag is an ^-dependent 
constant. In particular, the size of the pattern is Ag - 
times the size of the pattern for / = 1 with same A. The 
crossover length R* is expected to grows as y/N. 

For a background with £ > 1, the basis vectors at the 
unit cell are £e± and ££2, where ei, £2 are the basis vectors 
for £ — 1 background (see Eq. pf). Then the recipro- 
cal basis vectors are g\/£ and g2/£- From the observed 
patterns, we find that the line charge densities remain 
same for any £ (see Fig. 17 for an example of the patch 
boundaries). This implies that n\ and ri2 in eq. (8) 
are constrained to be multiples of I. Writing ri\ — Im, 
ri2 — In, we see that the patches can be labeled by the 
same pair of integers (m, n) as in the £ — 1 case, and the 



potential function (r) for general I is related to the 
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FIG. 17. An example of five patches meeting at a point, for 
a pattern on the background of Class I, £ = 2. It is easy to 
check that the line charge d ens ity for the vertical boundary 
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A = — 1/V3, same as in Fig. 

Fig. [TT] 



£ = 1 case by simple scaling: 



The color code is same as in 



W (r) = A e 4> {1) , (30) 

where is a scale factor. For ^ = 2, 3, 4 and 5 the 
values of Ag are approximately 2.34, 3.38, 4.41 and 5.37, 
respectively. We note that An increases approximately 
linearly with i. 



VI. NON-COMPACT PATTERNS WITH 
EXPONENT a < 1 

On the F-lattice, after some experimentation, we found 
that the background pattern having the periodicity of the 



tiling of plane with tilted rectangles, shown in Fig. 18 
produces patterns with interesting non-compact growth. 
We studied rectangles with aspect ratio I : (£ + 1), and the 
rectangles are tilted by 45° to the x-axis. Two such pe- 
riodic backgrounds are shown in Fig. [191 In these back- 



ground patterns, the sites with height zero, are arranged 
along the boundaries of tilted rectangles with two pos- 
sible orientations, and rest of the sites have height one. 
The stable height-patterns generated by adding N par- 
ticles and relaxing the configuration on these two back- 
grounds are shown in Fig. [20] and Fig. [21] respectively. 



The growing boundaries of the patches in the patterns 
are shown, in terms of the Q variables, in Fig. 22 and 
23 respectively. Again, we see that the patch boundaries 
are straight lines, with rational slopes. The plot of diam- 



eter 2A vs N, for these two patterns are shown in Fig. 24 
We see that the growth exponent a is approximately 0.6 
for figure 20 and 0.725 for figure 21 In general, value of 




FIG. 18. A schematic representation of the periodic tiling 
of the plane using tilted rectangles. Background height pat- 
terns with such periodicities on the F-lattice give rise to non- 
compact growth with the growth-exponent between 1/2 and 
1 



value 1 as density p Q of the background becomes close to 
1. 

There are unresolved areas of apparent solid color in 
the patterns, taking up a sizable fraction of the total area, 
e.g., two large regions of red color on both sides of Fig. 



22] In these regions, the pattern appears to be complex, 
suggesting either a large number of patch boundaries, or 
patches of non-zero areal excess charge density. However, 
the fractional area of these regions decreases with larger 
N. Also on comparing patterns with different /, we have 
seen that the fractional area of such regions decreases as I 
increases. A more detailed study of these patterns seems 
like an interesting problem for future investigations. 



VII. SUMMARY AND CONCLUDING 
REMARKS 

In this paper, we have studied two dimensional pat- 
terns formed in Abelian sandpile models by adding par- 
ticles at one site on an initial periodic background, where 
the diameter of the pattern grows as N a , with a > 1/2. 
Using some features observed in the pattern of adjacency 
of patches as an input, we are able to determine the exact 
asymptotic pattern in the specific case with a — 1, on a 
class I background. 

The patterns on class II backgrounds can also be char- 
acterized similarly. As noted earlier, some of the patches 
split into smaller parts. By using the 1/R transforma- 
tion, we can again determine the structure of the adja- 
cency graph. The graph for the pattern in Fig. ]8|b) is 
shown in Fig. 



25 



the exponent a is in range 1/2 < a < 1, and approaches 



It is a periodic lattice where half of 
the vertices of the hexagonal lattice are replaced by 3 ver- 
tices (colored in brown). The exact D- values for different 
patches can be easily determined. The determination of 
fm,n for this pattern then requires the solution of the 
Laplace's equation on this graph. It can be shown that 
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FIG. 19. The two backgrounds studied on the F-lattice. Unit 
cells of the periodic distribution of particles are shown by gray 
rectangular shades. The filled circles represent height and 
unfilled ones 1. 



a slight alteration of the graph, by drawing the missing 
edges in the small triangles shown in pink colors, does not 
change the pattern. Then the solution of the Laplace's 
equation can be reduced to the solution of a resistor net- 
work on this modified graph. The later can be further 
reduced to the resistor network on a hexagonal lattice, 
discussed by Atkinson et.al. [27], using the well-known 
Y — A transformation. We omit details of the analysis 
here. 

An important feature of the non-compact patterns is 
that, it can be characterized by a piece- wise linear func- 
tion. This characterization is simpler than that of the 
patterns with compact growth, where one requires piece- 
wise quadratic polynomials. We have shown that there 
are infinitely many backgrounds, on which the patterns 
have non-compact growth. It would be desirable to de- 
termine the exact value of a for different backgrounds 



Another interesting question is a possible connection of 
this problem to tropical algebra [29]. In tropical math- 
ematics, one defines operations similar to 'addition' and 
'multiplication' (denoted by © and © here) by 



a © b — max {a, b} 
a © b — a + 6, 



(31) 



where a, b are real parameters. Familiar properties of ad- 
dition and multiplication operators, like commutativity, 
associativity, existence of identity, distributive property 
continues to hold in the new definition. One can then 
define polynomials in several variables. The graph of a 
tropical polynomial is a piecewise linear function which is 
also convex. For example, consider the tropical function 



f(x) = a©x 2 ©6©x©c. 
In terms of standard algebra 

fix) — max {a + 2x, bx, c} 



(32) 



(33) 



The graph corresponding to this function is shown in Fig. 

m 

We note that for the pattern discussed in section 4, the 
potential function is piece-wise linear. It seems plausible 
that tropical polynomials may be useful to describe this 
function. In fact, tropical geometry have been discussed 
as possibly related to sandpile models [30, 3 1 J . For small 
values of N, our numerical study showed that is con- 
vex, if restricted to one sextant. However, for larger N, 
as shown in Fig. [27] we see that <fi is not convex even 
within one sextant. We conclude that it is not possible 
to represent the potential as a simple tropical polyno- 
mial. 
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Appendix: Relation to the theory of discrete 
analytic functions 

The sandpile patterns we studied are characterized in 
terms of discrete analytic functions (DAF) on different 
discretizations of the complex plane. For the pattern 
in Fig. 



11 it is the DAF on a hexagonal lattice, which 



increases logarithmically at large distances as in Eq. ( 27 ) 



showing non-compact growth studied in Sec. VI 



Studies of DAF started with the work of Kirchhoff on 
resistor networks [32Tf34~] . and has been studied subse- 
quently by many others [35] [36]. However, we have not 
encountered any work on DAF on many sheeted Riemann 
surfaces. In the following we present a way to determine 
DAF on a square discretization of Riemann surfaces. 
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FIG. 21. The pattern produced on the second background in Fig. 1 1 9 1 by adding N = 600 grains at a single site, and relaxing 
the configuration. Color code: White= 1 and Black= 0. Details can be viewed in the electronic version using zoom in. 
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FIG. 22. The pattern in terms of Q(r), showing the bound- 
aries of patches corresponding to Fig. [20] Color code: 
White= and Red=Non-zero. Details can be viewed in the 
electronic version using zoom in. 
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FIG. 24. The change in diameter as a function of N, for the 
patterns on the two backgrounds in Fig. |19| The numerical 
results fit well with straight lines of slope 0.6 and 0.725, for 
the backgrounds one and two, respectively. 




FIG. 23. The pattern in terms of Q(r), showing the bound- 
aries of patches corresponding to Fig. [2l] Color code: 
White= and Red=Non-zero. Details can be viewed in the 
electronic version using zoom in. 



Consider a square grid of points z — me + me, where 
m, n are integers and e is the lattice spacing. Let 
/ (me, ne) be a complex function defined at every site on 
the grid. The function / is defined to be discrete analytic 
[37] if it satisfies the discrete Cauchy Riemann condition 



f(z 3 )-f( Zl ) f{z A )-f{z 2 ) 



Z3 ~ Zi 



ZA - Z 2 



(a.i; 



at all elementary squares on the grid as shown in Fig. 28 
In complex analysis, simple examples of analytic func- 
tions are z n , and any polynomial of z n is also analytic 



For DAF, it is clear, using the linearity of equation (A.I ) 



that sum of DAF is also discrete analytic. However, not 
all positive integer powers of z are discrete analytic. It 



FIG. 25. The adjacency graph of the patches in the pattern 
in Fig. [8jb). The vertices corresponding to the brownish 
and greenish patches in the pattern (FigjTj) are denoted by 
different colors. 



is easy to check that the functions 1, 



are dis- 



crete analytic, but z 4 is not. We can however construct 
polynomial functions of Re(z) and Im(z), that are dis- 
crete analytic. Two such examples are z 4 — zze 2 and 
z 5 - 5z 2 ze 2 /2. 

We define a function F n (z,e) as a homogeneous poly- 
nomial in z, z and e, of degree n, which is a DAF. Using 
homogeneity, we have 



F n (z,e) = a n F n (- 



a a 



(A.2) 
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FIG. 26. Graph corresponding to the tropical function in 
equation (32 ) 



FIG. 28. A square grid on the complex plane. 



where 



- - - T> 




FIG. 27. (Color online.) Three dimension plot of the integer 
toppling function Tn (x, y) for a triangular pattern like in Fig. 
|11| but with N = 800. The plot shows a zoomed-in section 
in the region y > and y + \f?>x > 0. 



and then using a — e, we can express F n (z,e) in terms 
of F n (z, 1). This fixes F n (z,e) up to a multiplicative 
constant. The normalization is fixed by requiring that 
as e tends to zero, F n (z,e) — > z n . Then using Cauchy 
Riemann conditions it is easily seen that F n (z,e), for all 
integers n > 0, has a series expansion in e 2 of the form 



i 



n-3 V4 
7! 1 / n 



(4!)' 



10! 



,7 
1 / n 



(4!) 3 n- 9 \10 
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n — 7 



(A.4) 
(A.5) 

(A.6) 



and so on. For an integer n, this series will terminate 
after a finite number of terms, and all of them can be 
determined iteratively. 

It is possible to analytically continue the functions for 
rational values of n. For example, 



( n ) ( \ 



T(n + 1) 
4!r(n-2)' 



(A.7) 



Then, this analytic continuation of F n (z,e) provides us 
the discrete analytic functions which in the limit \z\ — > 
oo grows as z n , for any real positive values of n. It is 
interesting to note that the function D mjn , used in [23] to 
characterize the pattern in Fig. |2j^b) , is equal to Fi/ 2 (z — 
m + in, e — 1), up to a multiplicative constant. The 
patterns in the presence of a line of sinks, or near wedges 
studied in [2l] involve other rational values. For example, 
for the pattern near a line sink, one requires the function 
F 1/3 (z,l). 



F n (z,e) 



(A.3) 
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